$ontext
This code is for the following paper:
    Cai, Yongyang, William Brock, and Anastasios Xepapadeas (2022).
    Climate Change Impact on Economic Growth: Regional Climate Policy under Cooperation and Noncooperation.
    JAERE, forthcoming

Author: Yongyang Cai, The Ohio State University
Date: August 11, 2022
$offtext

$include Carbon_Calib.gms

set ii /North, Tropic, South/;

set thist(t) / 1900*2005 /;
set oo /Ocean/;

$offlisting

Table TA_RCP26(tt,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP26.csv
$offdelim;

Table TA_RCP45(tt,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP45.csv
$offdelim;

Table TA_RCP60(tt,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP60.csv
$offdelim;

Table TA_RCP85(tt,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP85.csv
$offdelim;

Table TLs(tt,n) 
$ondelim
$include OceanTempAnomaly.csv
$offdelim;

Table TL_hist(thist,oo) 
$ondelim
$include OceanTempAnomaly_historical.csv
$offdelim;

Table TA_hist(thist,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_historical.csv
$offdelim;

Table RFtotal(t,n) 
$ondelim
$include RCP_TotalRF.csv
$offdelim;

$onlisting
;

parameters 	
    fco22x   Forcings of equilibrium CO2 doubling (Wm-2)    /3.6813    /
    t2xco2   Equilibrium temp impact (oC per doubling CO2)    / 3.1  /
    
    xi1    
	xi2 
	xi3     
	B			

	xi4
    xi5
    xi7
    xi8
    
    TAs(ttt,n,ii)
    Forcings(ttt,n)
    Forcings_hist(thist)
;


Forcings(ttt,n) = RFtotal(ttt,n);
Forcings_hist(thist) = RFtotal(thist,'RCP85');

TAs(ttt,'RCP26',ii) = TA_RCP26(ttt,ii);
TAs(ttt,'RCP45',ii) = TA_RCP45(ttt,ii);
TAs(ttt,'RCP6',ii) = TA_RCP60(ttt,ii);
TAs(ttt,'RCP85',ii) = TA_RCP85(ttt,ii);
display Forcings, TAs;

****************************************************

positive variables
    epsTAT1(ttt,n,ii)
    epsTAT2(ttt,n,ii)
    epsTL1(ttt,n)
    epsTL2(ttt,n)
    epsTAT1hist(thist,ii)
    epsTAT2hist(thist,ii)
    epsTL1hist(thist)
    epsTL2hist(thist)    
;


VARIABLES
    xxi1    
	xxi2 
	xxi3     
	BB			

	xxi4
    xxi5
    xxi7
    xxi8
    
    TATM(ttt,n,ii)      Temperature of atmosphere in degrees C
    TOCEAN(ttt,n)       Temperature of lower oceans degrees C

    TATMhist(thist,ii)      Temperature of atmosphere in degrees C
    TOCEANhist(thist)       Temperature of lower oceans degrees C

    objTA
;

*****************************************************************

EQUATIONS

    TATMEQ_Tropic(ttt,n)        Temperature-climate equation for atmosphere
    TATMEQ_North(ttt,n)        Temperature-climate equation for atmosphere
    TATMEQ_South(ttt,n)
    TOCEANEQ(ttt,n)      Temperature-climate equation for lower oceans

    L1errTAT(ttt,n,ii)
    L1errTL(ttt,n)

    ClimSens
    Symmetric
    
    TATMEQ_Tropic_hist(thist)        Temperature-climate equation for atmosphere
    TATMEQ_North_hist(thist)        Temperature-climate equation for atmosphere
    TATMEQ_South_hist(thist)
    TOCEANEQ_hist(thist)      Temperature-climate equation for lower oceans

    L1errTAT_hist(thist,ii)
    L1errTL_hist(thist)
    
    L1objfunTA
;

** Equations of the model

TATMEQ_North(ttt+1,n)..   TATM(ttt+1,n,'North') =E= (1-BB)*TATM(ttt,n,'North') - 
						    xxi2*(TATM(ttt,n,'North')-TOCEAN(ttt,n)) -
						    xxi4*(TATM(ttt,n,'North')-TATM(ttt,n,'Tropic')) + 
						    (xxi1+xxi7)*Forcings(ttt,n);
TATMEQ_Tropic(ttt+1,n)..  TATM(ttt+1,n,'Tropic') =E= (1-BB)*TATM(ttt,n,'Tropic') - 
                            xxi2*(TATM(ttt,n,'Tropic')-TOCEAN(ttt,n)) - 
                            xxi4/2*(TATM(ttt,n,'Tropic')-TATM(ttt,n,'North')) - 
                            xxi5/2*(TATM(ttt,n,'Tropic')-TATM(ttt,n,'South')) +
                            (xxi1+xxi8)*Forcings(ttt,n);
TATMEQ_South(ttt+1,n)..   TATM(ttt+1,n,'South') =E= (1-BB)*TATM(ttt,n,'South') - 
						    xxi2*(TATM(ttt,n,'South')-TOCEAN(ttt,n)) -
						    xxi5*(TATM(ttt,n,'South')-TATM(ttt,n,'Tropic')) + 
						    xxi1*Forcings(ttt,n);

TOCEANEQ(ttt+1,n)..    TOCEAN(ttt+1,n) =E= TOCEAN(ttt,n) + sum(ii, xxi3*(TATM(ttt,n,ii)-TOCEAN(ttt,n))) +
                            xxi3*(TATM(ttt,n,'Tropic')-TOCEAN(ttt,n));

ClimSens..             BB*t2xco2 =e= (xxi1+(xxi7+2*xxi8)/4) * fco22x;
Symmetric..            xxi4 =e= xxi5;

L1errTAT(ttt,n,ii)..   TATM(ttt,n,ii)-TAs(ttt,n,ii) =E= epsTAT1(ttt,n,ii)-epsTAT2(ttt,n,ii);
L1errTL(ttt,n)..       TOCEAN(ttt,n)-TLs(ttt,n) =E= epsTL1(ttt,n)-epsTL2(ttt,n);

TATMEQ_Tropic_hist(thist+1)..  TATMhist(thist+1,'Tropic') =E= (1-BB)*TATMhist(thist,'Tropic') - 
                            xxi2*(TATMhist(thist,'Tropic')-TOCEANhist(thist)) - 
                            xxi4*(TATMhist(thist,'Tropic')-TATMhist(thist,'North')) - 
                            xxi5*(TATMhist(thist,'Tropic')-TATMhist(thist,'South')) +
                            (xxi1+xxi7)*Forcings_hist(thist);
TATMEQ_North_hist(thist+1)..   TATMhist(thist+1,'North') =E= (1-BB)*TATMhist(thist,'North') - 
						    xxi2*(TATMhist(thist,'North')-TOCEANhist(thist)) -
						    xxi4*(TATMhist(thist,'North')-TATMhist(thist,'Tropic')) + 
						    (xxi1+xxi8)*Forcings_hist(thist);
TATMEQ_South_hist(thist+1)..   TATMhist(thist+1,'South') =E= (1-BB)*TATMhist(thist,'South') - 
						    xxi2*(TATMhist(thist,'South')-TOCEANhist(thist)) -
						    xxi5*(TATMhist(thist,'South')-TATMhist(thist,'Tropic')) + 
						    xxi1*Forcings_hist(thist);

TOCEANEQ_hist(thist+1)..    TOCEANhist(thist+1) =E= TOCEANhist(thist) + sum(ii, xxi3*(TATMhist(thist,ii)-TOCEANhist(thist))) +
                            xxi3*(TATMhist(thist,'Tropic')-TOCEANhist(thist));


L1errTAT_hist(thist,ii)..   TATMhist(thist,ii)-TA_hist(thist,ii) =E= epsTAT1hist(thist,ii)-epsTAT2hist(thist,ii);
L1errTL_hist(thist)..       TOCEANhist(thist)-TL_hist(thist,'Ocean') =E= epsTL1hist(thist)-epsTL2hist(thist);


L1objfunTA..           objTA =e= SUM( (n,ttt,ii), (epsTAT1(ttt,n,ii)+epsTAT2(ttt,n,ii)) ) + 
				                 SUM((ttt,n), (epsTL1(ttt,n)+epsTL2(ttt,n))) +
				                 SUM( (thist,ii), (epsTAT1hist(thist,ii)+epsTAT2hist(thist,ii)) ) + 
				                 SUM(thist, (epsTL1hist(thist)+epsTL2hist(thist)));


**********************************

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = on;
option limrow = 0;
option limcol = 0;
option nlp = conopt;

*****************

xxi1.l = 0.1;
xxi2.l = 0.1;
xxi3.l = 0.002;
BB.l = 0.01;
xxi4.l = 0.4;
xxi5.l = 0.01;
xxi7.l = 0.01;
xxi8.l = 0.01;

xxi1.lo = 0.001;
xxi2.lo = 0.0005;
xxi3.lo = 0.0005;
xxi4.lo = 0;
xxi5.lo = 0.001;
xxi7.lo = 0;
xxi8.lo = 0;

BB.lo = 0.0001;
BB.up = 0.5;

TATM.l(ttt,n,ii)$(ord(ttt)=1) = TAs(ttt,n,ii);
TOCEAN.l(ttt,n)$(ord(ttt)=1) = TLs(ttt,n);
loop(ttt$(ord(ttt)<card(ttt)),
    TATM.l(ttt+1,n,'North') = (1-BB.l)*TATM.l(ttt,n,'North') - 
			xxi2.l*(TATM.l(ttt,n,'North')-TOCEAN.l(ttt,n)) - 
			xxi4.l*(TATM.l(ttt,n,'North')-TATM.l(ttt,n,'Tropic')) + 
			(xxi1.l+xxi7.l)*Forcings(ttt,n);
    TATM.l(ttt+1,n,'Tropic') = (1-BB.l)*TATM.l(ttt,n,'Tropic') - 
			xxi2.l*(TATM.l(ttt,n,'Tropic')-TOCEAN.l(ttt,n)) - 
			xxi4.l/2*(TATM.l(ttt,n,'Tropic')-TATM.l(ttt,n,'North')) - 
            xxi5.l/2*(TATM.l(ttt,n,'Tropic')-TATM.l(ttt,n,'South')) + 
			(xxi1.l+xxi8.l)*Forcings(ttt,n);
    TATM.l(ttt+1,n,'South') = (1-BB.l)*TATM.l(ttt,n,'South') - 
			xxi2.l*(TATM.l(ttt,n,'South')-TOCEAN.l(ttt,n)) - 
			xxi5.l*(TATM.l(ttt,n,'South')-TATM.l(ttt,n,'Tropic')) + 
			xxi1.l*Forcings(ttt,n);            
    TOCEAN.l(ttt+1,n) = TOCEAN.l(ttt,n) + sum(ii, xxi3.l*(TATM.l(ttt,n,ii)-TOCEAN.l(ttt,n))) +
            xxi3.l*(TATM.l(ttt,n,'Tropic')-TOCEAN.l(ttt,n));
);


TATM.fx(ttt,n,ii)$(ord(ttt)=1) = TAs(ttt,n,ii);
TOCEAN.fx(ttt,n)$(ord(ttt)=1) = TLs(ttt,n);
TATMhist.fx(thist,ii)$(ord(thist)=1) = 0;
TOCEANhist.fx(thist)$(ord(thist)=1) = 0;

model CalibTAL1 / TATMEQ_Tropic, TATMEQ_North, TATMEQ_South, TOCEANEQ, ClimSens, Symmetric,
    L1errTAT, L1errTL, TATMEQ_Tropic_hist, TATMEQ_North_hist, TATMEQ_South_hist, TOCEANEQ_hist,
    L1errTAT_hist, L1errTL_hist, L1objfunTA  /;

solve CalibTAL1 minimizing objTA using nlp ;

xi1 = xxi1.l;
xi2 = xxi2.l;
xi3 = xxi3.l;
B = BB.l;
xi4 = xxi4.l;
xi5 = xxi5.l;
xi7 = xxi7.l;
xi8 = xxi8.l

execute_unload 'climate_anomaly_param', b12, b23, xi1, xi2, xi3, B, xi4, xi5, xi7, xi8;

file SpatialAnomalyPaths;
put SpatialAnomalyPaths;
SpatialAnomalyPaths.nw = 12;
SpatialAnomalyPaths.nr = 2;
SpatialAnomalyPaths.nz = 1e-12;
loop((n,ttt),
  put ttt.tl:4:0;
  loop(ii, 
    put TATM.l(ttt,n,ii):14:6;
  );
  put TOCEAN.l(ttt,n):14:6;
  put /;
);
